fig.landings <- function(scenario){
	# Total landings, plot only if MPD loaded
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op <- par(no.readonly=T)
    SSplotCatch(opList[[1]][[4]],1)
    saveFig("fig.a")
    par(op)
  }else{
    cat(paste("\nError plotting 'Landings' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.b <- function(includeMPD=F,ylimit=6,useMaxYlim=T,opacity="20",scenario,...){
	# Spawning stock biomass

  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    mcbo <- A$mc$bo
    post.bo <- as.data.frame(window(mcmc(mcbo),start=Burn,thin=Thin))
    boci <- apply(post.bo,2,quantile,probs=c(0.025,0.5,0.975))

    mcbt <- A$mc.sbt
    post.bt <- as.data.frame(window(mcmc(mcbt),start=Burn,thin=Thin))
    btci <- apply(post.bt,2,quantile,probs=c(0.025,0.5,0.975))
    if(useMaxYlim){
      yUpperLimit <- max(btci,A$sbt)
    }else{
      yUpperLimit <- ylimit
    }
    matplot(A$yrs,t(btci),type="l",col=1,lty=c(2,1,2), lwd=2,ylim=c(0,yUpperLimit),ylab="Spawning biomass", las=1, ...)
    xx <- c(A$yrs,rev(A$yrs))
    yy <- c(btci[1,],rev(btci[3,]))
    shade <- getShade(1,opacity)
    polygon(xx,yy,density=NA,col=shade)
    points(A$yrs[1]-0.8,boci[2],col=1,pch=1)
    arrows(A$yrs[1]-0.8,boci[1],A$yrs[1]-0.8,boci[3],col=1, code=0, lwd=1.5)

    if(includeMPD && opList[[scenario]][[6]]){
      lines(A$yrs,A$sbt,type="l",col=mpdLineColor,lty=1, lwd=2,ylim=c(0,yUpperLimit), las=1, xlab="",ylab="")
      legend("topright",c("MPD estimate"),lty=1,lwd=2,col=mpdLineColor,bty="n")
    }
    saveFig("fig.b")
    par(op)
  }else{
    cat(paste("Error plotting 'Spawning Biomass' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.biomass.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
	# Spawning stock biomass
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],3)
    #if(useMaxYlim){
    #  yUpperLimit <- max(A$sbt)
    #}else{
    #  yUpperLimit <- ylimit
    #}
    #matplot(A$yrs,A$sbt,type="l",col=mpdLineColor,lty=1, lwd=2,ylim=c(0,yUpperLimit),ylab="Spawning biomass", las=1, ...)
    #points(A$yrs[1]-0.8,A$sbo,col=mpdLineColor,pch=1)
    saveFig("fig.biomassMPD")
    par(op)
  }else{
    cat(paste("Error plotting 'Spawning Biomass' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.acoustic.survey.index.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],2)
    saveFig("fig.acoustic.survey.index.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Acoustic Survey Index (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.acoustic.survey.log.index.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],5)
    saveFig("fig.acoustic.survey.log.index.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Acoustic Survey Log Index (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.acoustic.survey.index.oe.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],3)
    saveFig("fig.acoustic.survey.index.oe.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Acoustic Survey Index Observed vs. Expected (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.acoustic.survey.log.index.oe.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],6)
    saveFig("fig.acoustic.survey.log.index.oe.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Acoustic Survey Log Index Observed vs. Expected (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.acoustic.survey.log.index.oe.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotIndices(opList[[scenario]][[4]],6)
    saveFig("fig.acoustic.survey.log.index.oe.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Acoustic Survey Log Index Observed vs. Expected (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.end.year.expected.growth.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotBiology(opList[[scenario]][[4]],subplots=7)
    saveFig("fig.end.year.expected.growth.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Ending Year Expected Growth (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.log.recruitment.deviations.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotRecdevs(opList[[scenario]][[4]],2)
    saveFig("fig.log.recruitment.deviations.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Log Recruitment Deviations (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.recruitment.deviations.variance.check.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotRecdevs(opList[[scenario]][[4]],3)
    saveFig("fig.recruitment.deviations.variance.check.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Recruitment Deviations Variance Check (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.length.based.selectivities.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSelex(opList[[scenario]][[4]],subplot=1)
    saveFig("fig.length.based.selectivities.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Length Based Selectivities (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.age.based.selectivities.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSelex(opList[[scenario]][[4]],subplot=2)
    saveFig("fig.age.based.selectivities.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Age Based Selectivities (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.spawn.recruit.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSpawnrecruit(opList[[scenario]][[4]],subplot=1)
    saveFig("fig.spawn.recruit.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Spawn Recruit Curve (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.spawn.recruit.years.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSpawnrecruit(opList[[scenario]][[4]],subplot=2)
    saveFig("fig.spawn.recruit.years.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'Spawn Recruit Curve with Years (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.spr.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSPR(opList[[scenario]][[4]],subplots=1)
    saveFig("fig.spr.mpd")
    par(op)
  }else{
    cat(paste("Error plotting 'SPR (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.one.minus.spr.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSPR(opList[[scenario]][[4]],subplots=2)
    saveFig("fig.one.minus.spr.mpd")
    par(op)
  }else{
    cat(paste("Error plotting '1-SPR (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.one.minus.spr.over.spr40.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSPR(opList[[scenario]][[4]],subplots=3)
    saveFig("fig.one.minus.spr.over.spr40.mpd")
    par(op)
  }else{
    cat(paste("Error plotting '1-SPR/1-SPR40 (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.one.minus.spr.over.btarget.mpd <- function(ylimit=6,useMaxYlim=T,scenario,...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    SSplotSPR(opList[[scenario]][[4]],subplots=4)
    saveFig("fig.one.minus.spr.over.btarget.mpd")
    par(op)
  }else{
    cat(paste("Error plotting '1-SPR/1-SPR40 by B/Btarget (mpd)' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.c <- function(includeMPD=F,ylimit=3.5,useMaxYlim=T,opacity="20",scenario, ...){
	# Spawning stock depletion
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    mcdt <- A$mc.sbdepletion
    post.dt  <- as.data.frame(window(mcmc(mcdt),start=Burn,thin=Thin))
    dtci <- apply(post.dt,2,quantile,probs=c(0.025,0.5,0.975))
    if(useMaxYlim){
      yUpperLimit <- max(dtci,A$sbt/A$sbo)
    }else{
      yUpperLimit <- ylimit
    }
    matplot(A$yrs,t(dtci),type="l",col=1,lty=c(2,1,2), lwd=2,ylim=c(0,yUpperLimit),ylab="Spawning depletion", ...)
    xx <- c(A$yrs,rev(A$yrs))
    yy <- c(dtci[1,],rev(dtci[3,]))
    shade <- getShade(1,opacity)
    polygon(xx,yy,density=NA,col=shade)
    abline(h=0.40, lwd=mtLineWidth, col=mtLineColor, lty=mtLineType)
    if(includeMPD && opList[[scenario]][[6]]){
      lines(A$yrs,A$sbt/A$sbo,type="l",col=mpdLineColor,lty=1, lwd=2,ylim=c(0,yUpperLimit), las=1, xlab="",ylab="")
      legend("topright",c("Management target","MPD estimate"),lty=c(mtLineType,1),
             lwd=c(2,2),col=c(mtLineColor,mpdLineColor),bty="n")
    }else{
      legend("topright",c("Management target"),lty=c(mtLineType),lwd=c(mtLineWidth),col=c(mtLineColor),bty="n")
    }
    saveFig("fig.c")
    par(op)
  }else{
    cat(paste("\nError drawing 'Depletion' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

#this was figure e in 2010 assessment
fig.d <- function(ylimit=35,useMaxYlim=T, scenario,  ...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
	  #op	<- par(no.readonly=T)
    mc <- A$mc.rt
    mc.rt <- as.data.frame(window(mcmc(mc),start=Burn,thin=Thin))
    rt <- apply(mc.rt,2,quantile,probs=c(0.025,0.5,0.975)) #gets quantiles for number of age 1 recruits
    if(useMaxYlim){
      yUpperLimit <- max(rt)
    }else{
      yUpperLimit <- ylimit
    }
    xp <- plot(ryr, rt[2,], type="p", pch=20,ylim=c(0,yUpperLimit),ylab="Age-1 recruits (billions)", las=1, ...)
    arrows(ryr, rt[1, ],ryr,rt[3,],code=3,angle=90,length=0.01)
    ##points(xp,A$nt[,1],pch=19,cex=1)
    abline(h=median(as.matrix(mc.rt)),col=2)
    abline(h=mean(as.matrix(mc.rt)),col=3,lty=2)
    legend("topright",c("long-term median","long-term mean"),lty=c(1,2),pch=c(-1,-1),lwd=c(1,1),col=c(2,3),bty="n")
    saveFig("fig.d")
 	  #par(op)
  }else{
    cat(paste("\nError plotting 'Recruitment' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

#SPR status - fmsy
fig.e1 <- function(scenario){
	#The relative spawning potential ratio (1-spr)/(1-spr.at.msy)
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    spr <- A$mc.sprstatus_fmsy #read.table("ccam.spr",h=F)
    post.spr <- as.data.frame(window(mcmc(spr),start=Burn,thin=Thin))
    sprci <- apply(post.spr,2,quantile,probs=c(0.025,0.5,0.975))

    matplot(A$yr,t(sprci),type="l",col=c(2,1,2),lty=c(3,1,3), lwd=2, pch=c(-1, 0, 1),ylim=c(0,2)
            ,xlab="Year",ylab="(1-SPR)/(1-SPR at fmsy)")
    abline(h=1, lwd=mtLineWidth, col=mtLineColor, lty=mtLineType)
    text(1980, 1, "Management target", pos=3)

    saveFig("fige_fmsy")
    par(op)
  }else{
    cat(paste("\nError plotting 'SPRMSY' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

#SPR status - f40
fig.e2 <- function(scenario){
	#The relative spawning potential ratio (1-spr)/(1-spr.at.msy)
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    spr <- A$mc.sprstatus_f40 #read.table("ccam.spr",h=F)
    post.spr <- as.data.frame(window(mcmc(spr),start=Burn,thin=Thin))
    sprci <- apply(post.spr,2,quantile,probs=c(0.025,0.5,0.975))

    matplot(A$yr,t(sprci),type="l",col=c(2,1,2),lty=c(3,1,3), lwd=2, pch=c(-1, 0, 1),ylim=c(0,2)
            ,xlab="Year",ylab="(1-SPR)/(1-SPR at f40)")
    abline(h=1, lwd=mtLineWidth, col=mtLineColor, lty=mtLineType)
    text(1980, 1, "Management target", pos=3)
    saveFig("fige_f40")
    par(op)
  }else{
    cat(paste("\nError plotting 'SPRF40' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.g <- function(scenario){
	#The exploitation fraction (catch/3+biomass)
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    bt3 <- A$mc.bt3 #read.table("ccam.bt3",h=F)
    post.bt3 <- as.data.frame(window(mcmc(bt3),start=Burn,thin=Thin))
    bt3ci <- apply(post.bt3,2,quantile,probs=c(0.025,0.5,0.975))
    ct <- A$ct[1,]#/1.e6
    matplot(A$yr,ct/t(bt3ci[, 1:length(A$yr)]),type="l",col=c(2,1,2),lty=c(3,1,3), lwd=2, pch=c(-1, 0, 1)
            ,xlab="Year",ylab="Exploitation fraction (catch/3+ biomass)")

    saveFig("figg")
    par(op)
  }else{
    cat(paste("\nError plotting 'Exploitation fraction' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.h	<- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op <- par(no.readonly=T)
    par(mfcol=c(1,2))
    n <- length(A$yr)

    ft <- A$mc.ft #read.table("ccam.ft",h=F)
    sbt <- A$mc.sbt #read.table("ccam.sbt",h=F)
    mc <- A$mc #read.table("ccam.mcmc", h=T)
    sbstatus <- sbt/mc$bmsy
    fstatus <- ft/mc$fmsy

    sbci <- apply(sbstatus, 2, quantile, probs=c(0.5))
    ftci <- apply(fstatus, 2, quantile, probs=c(0.5))
    plot(sbci[1:n], ftci[1:n], type="n",xlim=c(0,1.2*max(sbci[1:n])), ylim=c(0,1.2*max(ftci)), xlab="Median Bt/Bmsy", ylab="Median Ft/Fmsy")
    abline(h=1,v=1,lty=2,col=2)
    lines(sbci[1:n], ftci[1:n],type="o")
    gletter(1)

	  #2nd plot
    sbstatus <- sbt/mc$bo
    sbci <- apply(sbstatus, 2, quantile, probs=c(0.5))
    ftci <- apply(ft, 2, quantile, probs=c(0.5))
    ssb <- seq(0,2,length=100)
    fp <- median(mc$fmsy)*(ssb-0.1)/0.3
    fp[ssb<=0.1] <- 0
    fp[ssb>0.4] <- median(mc$fmsy)
    maxY <- max(fp,A$ft)*1.1
    plot(ssb,fp,type="l",ylim=c(0,maxY),lwd=2,xlab="Median SBt/SBo",ylab="Median Ft")
    lines(ssb, fp, lwd=2)
    lines(sbci[1:n],ftci, type="o")
    gletter(2)

    saveFig("figh")
    par(op)
  }else{
    cat(paste("\nError plotting 'Bt/BmSY' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.i <- function(scenario){
	#plot the equilibrium yield curves
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op <- par(no.readonly=T)
    par(mfcol=c(2,2))
	  #A$equil comes from the ccam.rep file
    fe <- A$equil[, 2]
    ye <- A$equil[, 3]
    sde <- A$equil[, 5]
    spr <- A$equil[, 7]

    plot(fe, ye, type="l", xlab="Fishing mortality (Fe)", ylab="Equilibrium yield")
    gletter(1)

    plot(sde, ye, type="l", xlab="Spawning depletion", ylab="Equilibrium yield", lty=2, col=2)
    gletter(2)

    plot(spr,ye, type="l", xlab="Spawning potential ratio", ylab="Equilibrium yield", lty=3, col=3)
    gletter(3)

    matplot(cbind(fe, sde, spr), ye/max(ye)*100, type="l",xlab="Fe, depletion,  SPR",  ylab="Relative equilibrium yield")
    gletter(4)

    saveFig("figi")
    par(op)
  }else{
    cat(paste("\nError plotting 'Equilibrium yield' figure.  There are no MCMC output loaded for scenario",scenario,"\n\n"))
  }
}

fig.j <- function(scenario){
	#Relationship between fishing mortlaity ~ yield,  recruitment,  SBe,  SPR
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op <- par(no.readonly=T)
    par(mfcol=c(2,2))

    fe <- A$equil[, 2]
    ye <- A$equil[, 3]
    sde <- A$equil[, 5]
    re <- A$equil[, 6]
    spr <- A$equil[, 7]

    ix <- c(min(which(ye==max(ye))),min(which(sde<=0.4)) , min(which(spr<=0.4)))

    plot(fe, ye, type="l",xlab="", ylab="Equilibrium yield (million mt)", lwd=2)
    segments(fe[ix],0,fe[ix],ye[ix],lty=c(1, 2, 3))
    segments(0,ye[ix],fe[ix],ye[ix],lty=c(1, 2, 3))

    re <- re/re[1]
    plot(fe, re, type="l",xlab="", ylab="Relative recruitment", lwd=2)
    segments(fe[ix],0,fe[ix],re[ix],lty=c(1, 2, 3))
    segments(0,re[ix],fe[ix],re[ix],lty=c(1, 2, 3))

    plot(fe, sde, type="l",xlab="", ylab="Spawning depletion", lwd=2)
    segments(fe[ix],0,fe[ix],sde[ix],lty=c(1, 2, 3))
    segments(0,sde[ix],fe[ix],sde[ix],lty=c(1, 2, 3))

    plot(fe, spr, type="l",xlab="", ylab="Spawning Potential Ratio", ylim=c(0, 1), lwd=2)
    segments(fe[ix],0,fe[ix],spr[ix],lty=c(1, 2, 3))
    segments(0,spr[ix],fe[ix],spr[ix],lty=c(1, 2, 3))

    legend("topright", c("MSY", "SB40", "SPR40"), lty=1:3, bty="n")

    mtext("Equilibrium fishing mortality rate", 1, outer=T, line=-1)

    saveFig("figj")
    par(op)
  }else{
    cat(paste("\nError plotting 'Equilibrium F/Yield' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

## TABLES

table.b <- function() {
	mcbt <- A$mc.sbt
	post.bt <- as.data.frame(window(mcmc(mcbt),start=Burn,thin=Thin))
	btci <- t(apply(post.bt,2,quantile,probs=c(0.025,0.5,0.975)))

	filename <- paste(tabDir,"tableb_sbt.tex",sep="")
	filenamecsv <- paste(tabDir,"tableb_sbt.csv",sep="")
	t.b	<- cbind(A$yrs,btci[1:nyrs,])

	cap="Recent trends in estimated female spawning stock biomass (million mt) 
		based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","Female biomass")
	ncgrp	<- c(1,3)
	colnames(t.b)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.b, 10),file=filename,caption=cap,label="tableb",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.b,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.c <- function() {
	mcdt <- A$mc.sbdepletion
	post.dt  <- as.data.frame(window(mcmc(mcdt),start=Burn,thin=Thin))
	dtci <- t(apply(post.dt,2,quantile,probs=c(0.025,0.5,0.975)))

	filename<-paste(tabDir,"tablec_depletion.tex",sep="")
	filenamecsv<-paste(tabDir,"tablec_depletion.csv",sep="")
	t.c	<- cbind(A$yrs,dtci[1:nyrs,])

	cap="Recent trends in estimated spawning depletion level
		based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","Depletion")
	ncgrp	<- c(1,3)
	colnames(t.c)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.c, 10),file=filename,caption=cap,label="tablec",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.c,file=filenamecsv, row.names=F)
#  cat(paste("Saved table ",filename,"...\n",sep=""))
#  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.d <- function(){
	mcrt <- A$mc.rt
	post.rt  <- as.data.frame(window(mcmc(mcrt),start=Burn,thin=Thin))
	rtci <- t(apply(post.rt,2,quantile,probs=c(0.025,0.5,0.975)))

	filename<-paste(tabDir,"tabled_recruits.tex",sep="")
	filenamecsv<-paste(tabDir,"tabled_recruits.csv",sep="")
	t.d	<- cbind(ryr,rtci[1:length(ryr),])

	cap="Recent trends in estimated recruitment (billions of age 1 fish)
		based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","Recruits")
	ncgrp	<- c(1,3)
	colnames(t.d)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.d, 10),file=filename,caption=cap,label="tabled",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.d,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.e1 <- function() {
	spr <- A$mc.sprstatus_fmsy
	post.spr <- as.data.frame(window(mcmc(spr),start=Burn,thin=Thin))
	sprci <- t(apply(post.spr,2,quantile,probs=c(0.025,0.5,0.975)))

	filename<-paste(tabDir,"tablee_sprfmsy_status.tex",sep="")
	filenamecsv<-paste(tabDir,"tablee_sprfmsy_status.csv",sep="")
	t.e	<- cbind(A$yr,sprci)

	cap="Recent trends in (1-spr)/(1-spr at fmsy) based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","(1-spr)/(1-spr at fmsy)")
	ncgrp	<- c(1,3)
	colnames(t.e)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.e, 10),file=filename,caption=cap,label="tablee",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.e,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.e2 <- function(){
	spr <- A$mc.sprstatus_f40
	post.spr <- as.data.frame(window(mcmc(spr),start=Burn,thin=Thin))
	sprci <- t(apply(post.spr,2,quantile,probs=c(0.025,0.5,0.975)))

	filename<-paste(tabDir,"tablee_sprf40_status.tex",sep="")
	filenamecsv<-paste(tabDir,"tablee_sprf40_status.csv",sep="")
	t.e	<- cbind(A$yr,sprci)

	cap="Recent trends in (1-spr)/(1-spr at f40) based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","(1-spr)/(1-spr at f40)")
	ncgrp	<- c(1,3)
	colnames(t.e)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.e, 10),file=filename,caption=cap,label="tablee",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.e,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.f <- function() {
	bt3 <- A$mc.bt3
	ct <- A$ct[1,]
	cbt3 <- ct/bt3[, 1:length(A$yr)]
	post.cbt3 <- as.data.frame(window(mcmc(cbt3),start=Burn,thin=Thin))
	cbt3ci <- t(apply(post.cbt3,2,quantile,probs=c(0.025,0.5,0.975)))

	filename<-paste(tabDir,"tablef_exploitationFraction.tex",sep="")
	filenamecsv<-paste(tabDir,"tablef_exploitationFraction.csv",sep="")
	t.f	<- cbind(A$yr, cbt3ci)

	cap="Recent trends in Ct/Bt3+ based on 5000 systematic samples from the joint posterior distribution."
	cgrp <- c(" ","Ct/Bt3+")
	ncgrp	<- c(1,3)
	colnames(t.f)	<- c("Year",c("2.5%","median","97.5%"))

  dum <- latex(tail(t.f, 10),file=filename,caption=cap,label="tablef",rowname=NULL,cgroup=cgrp,n.cgroup=ncgrp)
  write.csv(t.f,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.h <- function(mle=T,tableType="ssb",perc=c(0.25,0.75),stock="Female",catchFactor=1e6,writeCSV=T){
  # tableDir:  absolute directory where the files will be stored
  # A:         an objects returned from reptolist, i.e. list of the rep file entries
  # Assumes A$mc.for exists and has the case sensetive column names: Year, CtStream, Rt, Sbt, f40spr, depletion, OY
  # streams: a vector of catch weights to include in the table. If a stream is less than zero, the OY will be used as the catch stream
  #          value in the table.  These f-based catch streams will be appended to the table in the order they were
  #          entered in the catchStreams vector.
  #          -They are found in the mcmc forecast file (tinss.for) in the CtStream column
  # mle: if True, use mle outputs.  If false, use MCMC output.
  # tableType = "ssb" means make table using spawning stock biomass (sbt column in A$mc.for) - makes es.table.h.1.csv
  # tableType = "depletion" means make table using reletive deletion (depletion column in A$mc.for) - makes es.table.h.2.csv
  # tableType = "f40spr" means make table using reletive spawning potential ratio (1-spr)/(1-0.4) (f40spr column in A$mc.for) - makes es.table.h.3.csv
  # stock is either "Female" or "All".  If Female the biomass is divided by 2.
  # catchFactor:  number to multiply catch by so it appears in the table in a nicer format
  # perc is a vector of the locations of the biomass vector to split the catch stream data on.

  #splits <- c(0,sort(perc),1)
  splits <- c(0,1)
  sections <- 1:(length(splits)-1)
  tmpMeds <- vector("numeric",length=length(sections))
  if(mle){
    forc <- A$mlefor
  }else{
    forc <- A$mcfor
  }
  years <- unique(forc$Year)
  streams <- unique(forc$CtStream)

  #medians <- as.data.frame(matrix(nrow=0,ncol=length(splits)+1))  # medians for the fixed catch streams
  #fmsyMedians <- as.data.frame(matrix(nrow=0,ncol=length(splits)+1)) # medians for the f-based catch streams
  #f40Medians <- as.data.frame(matrix(nrow=0,ncol=length(splits)+1)) # medians for the f-based catch streams
  #ssMedians <- as.data.frame(matrix(nrow=0,ncol=length(splits)+1)) # medians for the SS-OY-based catch streams

  medians <- as.data.frame(matrix(nrow=0,ncol=4))  # medians for the fixed catch streams
  fmsyMedians <- as.data.frame(matrix(nrow=0,ncol=4)) # medians for the f-based catch streams
  f40Medians <- as.data.frame(matrix(nrow=0,ncol=4)) # medians for the f-based catch streams
  ssMedians <- as.data.frame(matrix(nrow=0,ncol=4))# medians for the SS-OY-based catch streams

  for(year in years){
    for(stream in streams){
      yStream <- subset(forc,forc$Year==year & forc$CtStream==stream)
      specialStream <- median(as.numeric(yStream$OY))
      ctstream<-stream

      #colnames(medians) <- c("Year","CtStream",paste(splits[1:(length(splits)-1)]," - ",splits[2:length(splits)]))
      #colnames(fmsyMedians) <- c("Year","CtStream",paste(splits[1:(length(splits)-1)]," - ",splits[2:length(splits)]))
      #colnames(f40Medians) <- c("Year","CtStream",paste(splits[1:(length(splits)-1)]," - ",splits[2:length(splits)]))
      #colnames(ssMedians) <- c("Year","CtStream",paste(splits[1:(length(splits)-1)]," - ",splits[2:length(splits)]))

      colnames(medians) <- c("Year","CtStream","ABC","Median")
      colnames(fmsyMedians) <- c("Year","CtStream","ABC","Median")
      colnames(f40Medians) <- c("Year","CtStream","ABC","Median")
      colnames(ssMedians) <- c("Year","CtStream","ABC","Median")

      if(mle){
        sorted <- yStream
      }else{
        nsamp <- length(yStream[,1])
        yStream <- yStream[((Burn+1):nsamp),]
        sorted <- yStream #[order(yStream$Rt_2009),]
      }

      lenSorted <- nrow(sorted)
      for(section in sections){
        lims <- splits[section:(section+1)]
        lowerInd <- floor(lenSorted*as.numeric(lims[1]))
        upperInd <- floor(lenSorted*as.numeric(lims[2]))
        if(tableType=="ssb"){
          if(mle){
            tmpMeds[section] <- sorted$SBt
          }else{
            tmpMeds[section] <- median(as.numeric(sorted[lowerInd:upperInd,]$Sbt))
          }
        }
        if(tableType=="depletion") {
          if(mle){
            tmpMeds[section] <- sorted$depletion
          }else{
            tmpMeds[section] <- median(as.numeric(sorted[lowerInd:upperInd,]$depletion))
          }
        }
        if(tableType== "f40spr") {
	          if(mle){
	            tmpMeds[section] <- sorted$SPR40status
	          }else{
	            tmpMeds[section] <- median(as.numeric(sorted[lowerInd:upperInd,]$SPR40status))
	          }
        }
        if(section==length(sections)){
          if(stream>=0){
            medians <- rbind(medians,as.numeric(c(year,ctstream,as.numeric(stream)*catchFactor,tmpMeds)))
          }
          if(stream==.FMSYFORCFLAG){
            fmsyMedians <- rbind(fmsyMedians,as.numeric(c(year,ctstream,specialStream*catchFactor,tmpMeds)))
          }
          if(stream==.F40FORCFLAG){
            f40Medians <- rbind(f40Medians,as.numeric(c(year,ctstream,specialStream*catchFactor,tmpMeds)))
          }
         if(stream==.SSFORCFLAG){
           ssMedians <- rbind(ssMedians,c(year,ctstream,specialStream*catchFactor,tmpMeds))
          }
        }
      }
    }
  }

  sort.medians.for.output <- medians[order(medians$CtStream,medians$Year),]
  sort.medians.for.output <- rbind(sort.medians.for.output,fmsyMedians)  # add fmsy-based catch streams
  sort.medians.for.output <- rbind(sort.medians.for.output,f40Medians)  # add f40-based catch streams
  sort.medians.for.output <- rbind(sort.medians.for.output,ssMedians)  # add SS-OY-based catch streams

  if(writeCSV){
    if(tableType=="ssb")
      fn <- paste(tabDir,"table.h1.ssb",sep="")
    if(tableType=="depletion")
      fn <- paste(tabDir,"table.h2.depletion.",sep="")
    if(tableType=="f40spr")
      fn <- paste(tabDir,"table.h3.f40spr",sep="")
    if(mle){
      fn <- paste(fn,"_mle.csv",sep="")
    }else{
      fn <- paste(fn,"_mcmc.csv",sep="")
    }
    # remove NAs from output
    sort.medians.for.output <- sort.medians.for.output[!is.na(sort.medians.for.output$CtStream),]
    sort.medians.for.output <- sort.medians.for.output[!is.na(sort.medians.for.output$Year),]
    write.csv(sort.medians.for.output,file=fn,row.names=F)
    cat(paste("Saved table ",fn,"...\n",sep=""))
  }
}

table.i <- function(){
	#This is the final summary table for the executive summary
	lbl <- c("Unfished SBo (million mt)", 
           "Unfished age-1 recruits (billions)", 
           "\\underline{\\emph{\\textbf{REFERENCE POINTS based on SB$_{40\\%}$}}}", 
           "MSY proxy spawning biomass SB$_{40\\%}$", 
           "SPR resulting in SB$_{40\\%}$", 
           "Exploitation fraction (ct/Bt3) resulting in SB$_{40\\%}$", 
           "Yield with SB$_{40\\%}$", 
           "\\underline{\\emph{\\textbf{REFERENCE POINTS based on SPR$_{40\\%}$}}}", 
           "Spawning biomass at SPR$_{40\\%}$", 
           "SPR", 
           "Exploitation fraction (ct/Bt3) resulting in SPR$_{40\\%}$", 
           "Yield with SPR$_{40\\%}$", 
           "\\underline{\\emph{\\textbf{REFERENCE POINTS based on MSY}}}", 
           "Spawning biomass at MSY", 
           "SPR at MSY", 
           "Exploitation fraction (ct/Bt3) at MSY", 
           "MSY")

  mcs <- A$mcRefPoints
  post.mcs <- as.data.frame(window(mcmc(mcs),start=Burn,thin=Thin))
  mcsci <- t(apply(post.mcs,2,quantile,probs=c(0.025,0.5,0.975)))
  ti <- mcsci
  ti <- rbind(ti[1:2, ], rep("", 3), ti[3:6, ], rep("", 3), ti[8:11, ], rep("", 3), ti[12:15, ])
  ti <- cbind(lbl, (ti))
  colnames(ti) <- c("Quantity", "2.5\\% percentile", "Median", "97.5\\% percentile")
  fn <- paste(tabDir,"table.i.csv",sep="")
  write.csv(ti,file=fn, row.names=F)
}

table.g1 <- function(){
  s <- read.csv(file="HakeManagement.csv")
  cap <- "Recent trend in Hake Management performance."
	filename <- paste(tabDir,"tablec.tex",sep="")
  dimnames(s)[[2]] <- c("Year","Landings","OY(mt)","ABC(mt)","Landings/OY(\\%)")
	latex(s,file=filename,caption=cap,label="tablec",rowname=NULL)
}

table.g2<-function(){
	#This is a sideways table
  s <- read.csv(file="Summarytable.CSV")
	cap <- "Summary of recent trends in Pacific hake exploitation and stock levels."
	filename <- paste(tabDir,"tableg.tex",sep="")
  latex(s[, c(1, 2:11)],file=filename,caption=cap,label="tableg",rowname=NULL)
}
